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We propose a nonperturbative quantum dissipation theory, in term of hierarchical quantum mas- 
ter equation. It may be used with a great degree of confidence to various dynamics systems in 
condensed phases. The theoretical development is rooted in an improved semiclassical treatment of 
Drude bath, beyond the conventional high temperature approximations. It leads to the new theory a 
simple modification but important improvement over the conventional stochastic Liouville equation 
theory, without extra numerical cost. Its broad range of validity and applicability is extensively 
demonstrated with two-level electron transfer model systems, where the new theory can be con- 
sidered as the modified Zusman equation. We also present a criterion, which depends only on the 
system-bath coupling strength, characteristic bath memory time, and temperature, to estimate the 
performance of the hierarchical quantum master equation. 



I. INTRODUCTION 

Dissipation is often inevitable and plays important 
roles in many systems. The key quantity in quantum dis- 
sipation theory (QDT) is the reduced system density op- 
erator, p{t) = trBPtot(i); i-e., the trace of total composite 
one over bath subspace. For Gaussian bath, exact QDT 
can be formulated with path integralii^i^ or its differen- 
tial version in terms of hierarchical equations of motion 
(HEOM)4>^>&i>a>aiioaiii2^^a5a6^ However, exact ap- 
proaches are numerically expensive in general. In this 
paper, we propose an approximate HEOM theory, which 
will be termed hereafter as hierarchical quantum master 
equation (HQME) . Compared with the traditional high- 
temperature approximation (HTA) schemes such as the 
stochastic Liouville equation^ '^^i^^'^°i^^ the new method 
requires about the same numerical effort, but enjoys a 
greatly improved range of validity. Moreover, the HQME 
supports also a convenient and versatile criterion of ap- 
plicability that seems to be rather insensitive to specific 
systems. 

We will exemplify the validity and applicability of the 
present theory with the standard electron transfer spin- 
boson model system. In this case, the HQME is ana- 
lytically solvable, and its high-temperature limit is just 
the celebrated Zusman equation The ZE 

treats the effect of bath via a diffusive solvation co- 
ordinate. Its validity has been a subject of study for 
years»^'29^>^^^^ With the aid of analytical so- 
lutions, we can now readily exploit the range of validity 
for both HQME and ZE, over the full parameters space 
of electron transfer systems. 

The present development is rooted in an improved 
semiclassical treatment of the fluctuation-dissipation 
theorem. Consider a Gaussian stochastic bath variable 
Fb(<). It can be the fluctuating solvation coordinate 
U{t) — {U)b in electron transfer systems, or the fluctuat- 



ing transition frequency 5L0eg{t)— {Sujeg) ^ in spectroscopy. 
Its effect on system is completely described by the bath 
correlation function, C{t) = {F-B{t)Fs{0))-B- The classical 
Gaussian-Markovian description assumes 



(1) 



Adopted is the classical fluctuation-dissipation relation, 
such as {U'^)b-{U)1 w 2kBT{U)B ^ iXksT. This model 
has been widely used, for example, in the spectroscopic 
motional narrowing problems J^ii^i^ It also leads to the 
stochastic Liouville equation description of reduced sys- 
tem dynamics However, it discards the imagi- 
nary part that is responsible for the spectroscopic Stokes 
shift or solvent reorganization. In other words, the clas- 
sical bath correlation function does not consider the back 
action of system on bath. The conventional HTA scheme 
adopts 



CHTA(t) = A(2fcsr-z7)e-'^*. 



(2) 



Here, the imaginary part assumes no approximation in 
the Drude dissipation model, while the real part remains 
its classical form. This scheme does account for the back 
action of system on bath. However, the reduced system 
density matrix dynamics based on CHTA(i) encounters 
rather often the positivity violation problem that was 
originally not suffered by that based on classical Cc\{t). 

The above two conventional schemes are the low-order 
approximations of fluctuation-dissipation theorem, to- 
gether with the Drude bath spectral density model. 



up -I- 7^ 

The exact fluctuation-dissipation theorem reads^ 

J(u;) 



C(t) = - / dtje-"^*- 

TI" J-oo 1 



(3) 



(4) 
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with j3 = l/lksT) being the inverse temperature. Con- 
sider the bosonic function in the expansion: 

The classical Cd (t) [Eq. ([1])] uses only the lowest order ex- 
pansion; the high-temperature Chta(^) [Eq. ([2])] includes 
also the second term. 

The proposed HQME approach is based on the follow- 
ing semiclassical bath correlation function, 

Cscit) - CHx.(t) - ^[7e-^* - 26it)]. (6) 

It uses all the three lowest terms of Eq. ([S]), thus, im- 
proving over CHTA(i) by two orders in Plu. Interestingly, 
the resulting HQME for reduced dynamics costs almost 
no additional numerical effort, but largely overcomes the 
aforementioned positivity problem. This remarkable fea- 
ture will be exemplified in simple electron transfer sys- 
tems, where the HTA limit of HQME has been shown to 
be equivalent to the ZE.— 

The remainder of paper is organized as follows. We 
present the HQME on the basis of Csc {t) [Eq. (O] , its 
stochastic Liouville equation description, and its contin- 
ued fraction Green's function theory, respectively, in the 
three subsections of Sec. [Hi The exact HEOM formal- 
ism based on C (t) of Eq. ^ is briefed in Appendix. In 
Sec. mil we consider the two-level electron transfer spin- 
boson model, in which the continued fraction Green's 
function theory of HQME can be analytically resolved 
and also the ZE is recovered in the HTA limit. Ex- 
emplified with this model, numerical studies on validity 
and applicability of HQME are carried out in Sec. lIVI 
We show that the HQME remarkably improves over its 
HTA/ZE counterpart, in terms of both positivity and 
electron transfer dynamics. Comments and discussions 
about the criterion for the applicability of HQME to ar- 
bitrary systems are presented in Sec.|Vl Finally, wc con- 
clude the paper in Sec.lVIl 



where 

Cr^2XkBT~-fA, c, = A7, A=-^. (9) 

bKB-L 

The corresponding HQME can then be constructed via 
the standard calculus-on-path-integral algebra i^i^°'^-^ It 
reads 

Pn{t) = -{iC + 5TZ + n'y)pn{t) 

- i^/nApn-i{t) - iVn + lBpn+i{t), (10) 

with 

6n6 = A[Q,[Q,d]], 6d = /M[Q,0], (lla) 
AO^ {cr[Q,d]-ic,{Q,d})/y^\. (lib) 

In this formalism, the reduced system density operator, 
p{t) = po{t), couples hierarchically with a set of well- 
defined auxiliary density operators {pn>o(0}- The hier- 
archical construction resolves not just system-bath cou- 
pling strengths but also memory time scales. Each p„, 
having been scaled by the factor (n!|cr|")~^/^, if com- 
pared with Ref.[l3 or Ref.lH, is now dimensionless and 
possesses a unified error tolerance as that of po- Thus, 
an efficient on-the-fly filtering algorithm that also auto- 
matically truncates the hierarchy can be applied.^'* 

Note that the exact HEOM formalism can be 
constructed in principle for arbitrary non-Markovian 
dissipation;^^'^- see Appendix for the exact theory of 
Drude dissipation. There the bath correlation function 
C{t) is expanded in Matsubara exponential series of the 
exact fiuctuation-dissipation theorem [Eq. ^] . The ex- 
act theory is generally expensive, even with the state-of- 
the-art numerical filtering algorithmji^ Apparently, the 
HQME [Eq. pH)) ] is numerically appealing, much more 
practical to large systems, if the effect of involving ap- 
proximation can be assessed in advance. We shall discuss 
this issue later in Sec.lIVI and Sec.lVl 



II. THEORY 

A. Hierarchical quantum master equation 

Consider the total Hamiltonian in the form of 

H^(t)=H + QF^{t). (7) 

Denote CO = [iJ, O] as the reduced system Liouvillian, 
and set Ti = 1 throughout this paper. The last term of 
Eq. ([7]) is the system-bath coupling, in which the system 
operator Q defines the dissipative mode, through which 
the stochastic bath operator or generalized Langevin 
force Fs{t) acting on the system. Recast the semiclas- 
sical bath correlation function of Eq. ^ as 

Csc{t) = {cr - ic^)e'^' + 2M{t), (8) 



B. Stochastic Liouville equation description 

Note that the HQME [Eq. HU])] has basically the 
same mathematical form, and thus about the same nu- 
merical cost, as the conventional stochastic Liouville 
equation.— The latter is based on the classical 
bath Cc\{t) [Eq. ([I])], and can be recovered from Eq. (fTU]) 
by setting Ci = A = 0. Thus, the HQME supports 
the same physical picture, as described by t) in the 
stochastic Liouville equation, with the diffusive solva- 
tion variable O being introduced for the effect of bath. 
The reduced system density operator is evaluated via 
p{t) = / dnp{n,t). The HQME [Eq. ([TOl)] resolves the 
stochastic description as 

oo 

p(r!,t) = ^p„(t)0„(i7), (12) 

n=0 
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where (/)„(0) e" — ^'^''(17), with ^^'^''(^2) being 

the normaHzed harmonic eigenfunction, is the right- 
eigenfunction of the diffusion operator;^ 



d_ 

on 



(13) 



With the same algebra of Ref.lla, we can show that the 
HQME [Eq. (fTU]) ] is equivalent to the following stochastic 
Liouville equation description, 



di 




In the classical bath limit (setting A = q = 0), the 
above equation reduces to the conventional stochastic Li- 
ouville equation4iii^ii2ii^i2i We have therefore extended 
the stochastic Liouville equation to not just the HTA 
bath, but also the present improved semiclassical scheme. 
Moreover, for the spin-boson system, the HTA (A = 0) 
version of Eq. has been recently showi*i^ to be iden- 
tical to the Zusman equation (ZE) i^i^^i^iiS^ Therefore, 
Eq. pop or Eq. can also be considered as a general- 
ized and modified ZE to arbitrary systems, with much 
improved validity range of parameters; see Sec. lIVI for a 
thorough demonstration. 



C. Continued fraction Green's function formalism 

The HQME formalism can in general apply to the sys- 
tems in the presence of external time-dependent field 
driving. In this case, the initial conditions to Eq. (fTU)) 
are the steady-states {/?„(< = 0) = p^} before the exter- 
nal field driving. This initial conditions can be evaluated 
by setting p„ = 0, leading to Eq. (fTO| a set of linear 
equations, under the constraint of Trpo = 1. It results in 
Pn>o{i = 0) 7^ generally, due to the initial system-bath 
coupling. 

For the population transfer systems to be studied in 
the absence of external field driving, we set the initial 
state to be pn{t = 0) = p(0)(5„o. This corresponds to the 
initial total density matrix factorization ansatz. In this 
case, the HQME [Eq. (fTU]) ] can be formally resolved with 
a continued fraction Green's function formalism4^'^ '^^i'^^ 
Following the same procedure as Ref. HI, we introduce 
the hierarchical Liouville-space propagators {lAn>o{t)} 



p„{t) = e-"^*Z^„(t)p(0); with Un{0) = <5„o, (15) 



and recast Eq. (fTO]) as (setting iC = iC + STZ) 

Unit) = -iC'Un{t) - iVnAe''*Un-i{t) 

-iy/n+lBe-^^Un+i{t). (16) 
In the Laplace domain, we have 

<5„o = (s + iC')Un{s) + iy/nJ\Un-i{s - 7) 

+ zV^mSZ^„+i(s + 7)- (17) 

Define the hierarchical Liouville-space Green's functions 
{e(")(s)} via 

U^{s) = g^''\s), (18a) 
W„(s) = -i\AI^(")(s)^„-i(s-7); n>0. (18b) 
Then Eq. dTT]) leads to 

1 



with 



H(") (s) = {n + l)ee("+i) (s + -f)A . 



(19a) 



(19b) 



The above equations constitute the continued fraction 
formalism to evaluate each individual H("^(s) or ^/''"-'(s). 

Note that Q''°\s) = Q[s) and its associated H(°)(s) = 
n(s) are the primary Green's function and dissipation 
kernel resolution, respectively. The reduced density op- 
erator p{s) = Q{s)p{0) satisfies 

sp{s) - p(0) - -iC'p{s) - U(s)p{s), (20) 

which in the time domain reads 



m = -ic'pit) 



dTfl{t-T)p{T). 



(21) 



This is nonperturbative quantum master equation, with 
II(< — t) [or H(s)] evaluated nonperturbatively via the 
continued fraction formalism and iC = iC + STZ. 

To obtain the kinetic rate equations, we start with 
Eq. (|20p and formally eliminate the coherence compo- 
nents of the reduced density matrix, resulting in 



sP{s)-P{0) = K{s)P{s). 
Its time-domain counterpart reads 



I. Jo 



(22) 



(23) 



where Kjk{s) are transfer rates resolutions for the transi- 
tion from the state k to the state j. They can be obtained 
via 



(24) 



Here, Tpc, Tcci Tcp, and Tpp denote the coherence- 
to-population, coherence-to-coherence, population-to- 
coherence, and population-to-population transfer matri- 
ces, respectively, defined by Eq. (j20p in a given represen- 
tation. 
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III. ANALYTICAL RESOLUTION OF 
ELECTRON TRANSFER DYNAMICS 

Consider hereafter the standard electron transfer 
model system. The total system-plus-bath composite 
Hamiltonian assumes Ht = ha\a){a\ + {hb + E°)\b){b\ + 
V{\a){b\ + \b){a\). Here, E° denotes the reaction en- 
dothermicity, and V is the transfer coupling matrix ele- 
ment that is assumed independent of the solvent degrees 
of freedom; ha or hh is the solvent Hamiltonian for the 
system in the donor \a) or acceptor \b) state, respectively. 
Their difference defines the solvation coordinate. 



U ^ hb- ha- 



(25) 



The system is assumed to be initially in the donor state, 
with pt(0) — \a){a\p'^, where the bath Hamiltonian as- 
sumes hs — ha, i.e., p'S^ = Pa'^ oc e~^'*". ft also de- 
fines the solvation reorganization energy A — {U)b = 
truiUp'^) and the generalized Langevin force FB{t) = 
gth-Bt ^jj _ ^-^g-ih-Bt physically corresponds to the dif- 
fusive variable fl introduced in Sec. lHBI 

The reduced electron transfer system Hamiltonian 
reads now 

H={E° + X}\b){b\ + V{\a){b\ + \b){a\), (26) 

while the dissipative system mode is 

Q^\b){b\. (27) 

The following derivation of the analytical solutions to the 
transfer dynamics follows the same algebra as Ref.lSGl. 
where the HTA version was treated. By analyzing the 
tensor elements involved in Eq. (|19bp for the specified 
dissipative mode Q of Eq. l|77|) , we find that the only 
nonzero tensor elements of H^"^ remain to be 



H 



ba,ba^ 



y 



in) 



t(") Jn) 



H 



(n) 
ba,bb^ 



(28) 



and their Hermitian conjugate counterparts. They are 
related to the Green's function tensor elements, 

Yin) ^ Qin) yin) ^ Qin) ^(„) ^ ^(n) 
^ — ^ba,ba' ~ ^ba.ab^ ^ — ^ba,bb' ^^-^i 

via [cf. Eq. (|19bp and denoting rj = Cr — ici] 

x(")(s) = ?7(n-f l)X("+i)(s-H7), (30a) 
2/(")(s) = -r7*(n+l)r("+i)(s-K7), (30b) 
z(")(s) = (77-ry*)(n + l)z("+i)(s + 7). (30c) 



From Eq. (|19ap and the associated Dyson equation 
method)^ we then obtain 

X(")(s) = [a(")(s) + /3(")(5)]7C^"Hs), (31a) 
r(")(s) - [/5(")(s) - y^"\s)]/<:^"Hs), (31b) 

Z(")(s) = -i([z(")(s) - iV]X^''Hs) 
s I 

-f [z(")(s)-il/]*y(")(s)}> (31c) 



with (^"Hs) 
and 



+ /3(»)(s)|2 _ |/?(")(s) - y(")(s)|2 



a("Hs) = s + i{E° + A) + A -I- a;(")(s), (32a) 
/3(")(s) EE s-^V[2V + iz(")(s)]. (32b) 

The above formulations [Eqs. ((28|) - ((32)) ] constitute the in- 
verse recursive analytical evaluation of the required dissi- 
pation kernel H(s) = H'^''^(s), and also the Green's func- 
tion. 

The forward and backward rate resolutions k{s) and 
k'{s) for the present two-level system of study can then 
be evaluated by using Eq. ([24| , resulting in 



k{s) = 2T/2Re 



a{s) +v{s) 
\a{sW-\y{sW 



and 



k'{s) = 2V'^Re 



[a[s)+y{s)][l-iz*{s)lV] 
K.)P-|y(,s)|2 



(33a) 



(33b) 



The rate constants are the values at s = for their steady 
state nature. The equilibrium reduced density matrix is 
obtained by solving [iC -I- H(s = 0)]p'^'^ = 0, together 
with Trp'"! = 1. The solutions with s = are 



Pbb ^ ^ Paa — 



cq cq 
Pab = Pba = 



Re (a -I- y) 



Re[{a + y){2-iz*/V)y 
Rc z 



Rc[{a + y){2-iz*/V)y 



(34a) 
(34b) 



Note that the electron transfer system considered here 
is different from the standard spin-boson model by their 
initial equilibrium bath states. The former is determined 
by hs = ha, while the latter by hg = ■^{ha + hi,). This 
distinction results in different forms of the reduced sys- 
tem Hamiltonian and dissipative mode; thus different 
p(t) dynamics. The steady-state behaviors would be the 
same, as required by thermodynamics principles, if the 
exact QDT is used. However, as here approximations 
are involved in treating the bath correlation function, 
the steady state behaviors such as rate constants and 
equilibrium system density matrix can be different in the 
aforementioned two model systems. Nevertheless, for ei- 
ther the standard electron transfer or spin-boson system, 
the algebra to analytical solutions is same. 



IV. NUMERICAL VALIDATIONS 

In this section we shall show that the proposed HQME 
is remarkably superior over the original HTA/ZE scheme. 
We address the issues of validity and applicability by con- 
sidering the positivity and accuracy of the reduced den- 
sity matrix dynamics. A versatile criterion for the range 
of applicability of the HQME will be constructed later; 
see the next section. Numerical results are all reported 
in unit of fcsT' = 1. 
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FIG. 1: Positivity diagram of Zusman equation, where the 
P-region is indicated by the arrow or arrows associated with 
each curve or each pair of curves in same color, (a): Dia- 
gram in (A, 7; V^)-subspace (with some selected values of V) 
for classical barrierless {E° + A = 0) and symmetric {E° — 0; 
inset) systems, (b): Diagram in (A, 7; i?°)-subspace for V — 1 
and 0.1 (inset), (c): Same as (b) but plotted in (i5°,A;7)- 
subspace. Unit of fcsT = 1 is used. 



For the issue of positivity, we focus on the asymptotic 
regime, characterized by the so-called P-region in the 
electron transfer parameters space, where the rate con- 
stants, k = fc(s — 0) and k = k'{s — 0), and the equi- 
librium density matrix all satisfy the positivity require- 
ment. The latter amounts to p'^pl'^ > p'^pla- With 
the aid of the analytic results presented in Eqs. ([55]) 
and (j34p . we explore thoroughly the P-regions in the 
temperature-scaled parameters ( A, 7, , y)-space of the 
electron transfer system. 

The resulting positivity diagrams of the ZE are re- 
ported in Fig.[TJ identical for endothermicity E° and 
—E°. The P-region is indicated by the arrow(s) associ- 
ated with each curve (or pair of curves of same color). In 
general, the P-region is larger for higher temperature, as 
expected. Also, the P-region itself does not vary when 
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FIG. 2: Population evolutions: Exact (black), HQME (red), 
and HTA/ZE (blue), for the specified systems. For the pa- 
rameter K, see Eq. (|35|l . Unit of ksT = 1 is used. 



=0.95 
Y=7 




200 



FIG. 3: Same as Fig. (2] but for systems with the common 
{E° , V, A) = (—0.5, 1, 13), and different values of 7, specified 
individually in each panel. Unit of fcsT = 1 is used. 



the transfer coupling strength (3V < 0.1, but changes 
quite dramatically as (3V increases. Observations for two 
specific systems, as depicted in Fig.[TIa) and its inset, are 
as follows. For a classical barrierless {E° + A = 0) sys- 
tem, the P-region in (A, 7)-subspace quickly shrinks with 
PV > 0.1, and the ZE violates positivity completely for 
f3V > 1.04. For a symmetric {E° = 0) system, the P- 
region does not change much, if (3V < 1, covering over all 
/?A < 11; however, when 1 < (3V < 4.5, it confines only 
within certain range of (3X depending on the value of f3"/. 
The ZE violates positivity completely for V > 4.5 for a 



6 



symmetric system. Apparently, the applicability range 
of HTA/ZE scheme depends not just on bath interaction 
parameters, but also sensitively on system itself, lead- 
ing to the complicated P-region diagrams, as depicted in 
Fig.ffl 

The HQME is simply superb. While the ZE is sub- 
ject to severe positivity violation especially for PV > 1, 
the new scheme is found to have P-region covering over 
a broad range of parameters space, as tested: (3X < 25, 
/37 < 100, \I3E°\ < 25, and PV < 100. In other words, 
the HQME scheme preserves positivity, at least asymp- 
totically, in almost entire parameters space of practical 
interest. Moreover, it is likely to support a convenient 
and system-insensitive criterion for its range of applica- 
bility; see the next section. 

The HQME shows also its superiority in time evolu- 
tion. It follows the exact results closer than the HTA/ZE 
does. Figure [2] depicts the transfer population evolutions 
of HQME, ZE, and exact HEOM. Chosen here are a clas- 
sical barrierless system {E° + A = 0) in the main panel 
and a symmetric system {E° = 0) in the inset. In con- 
trast to the ZE that evolves into unphysical regime of 
positivity violation, the HQME remains accurate quan- 
titatively and semi-quantitatively, respectively, for the 
two specified systems in study. Figure [3] is for three 
other systems, with common (3V, (3X, and f3E°, but dif- 
ferent values of /3j, where the intermediate one (/37 = 5) 
leads to the ZE dynamics violation of positivity. In all 
cases studied, the HQME performs much better than the 
HTA/ZE does. The superiority can be qualitative. While 
the HTA/ZE violates the positivity, the HQME can re- 
main even quantitatively applicable. 



V. DISCUSSION AND COMMENTS 
A. Criterion and measure on applicability 

For the issue of applicability, we propose to use 

K = V6r(7)/(/3A7), (35a) 

with 

/3r(7) - 6 + v/12 + (/37)' , (35b) 

to estimate both the range of applicability (setting to be 
K > 1) and the quality of HQME dynamics. The HQME 
dynamics tends to be more accurate for a larger k case. 
Justifications on the above criterion will be made in the 
next subsection. 

We had indicated in Fig. [2] and Fig. [3] their associ- 
ated values of the k parameter. The above system- 
independent criterion based on n is shown to be overall 
satisfactory. A quantitative agreement of HQME with 
the exact HEOM dynamics would depend on the specific 
details of systems. However, the system dependence is 
relatively insensitive, comparing to that on the nature of 
bath characterized by the parameter k. 
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FIG. 4: Population evolutions; Exact (black), HQME (red), 
and M-HQME (green), for {E°,V) = (-1,0.5) (lower) and 
{E° , V) = (—5, 3.5) (upper), respectively, with same bath pa- 
rameters of (A, 7) = (20, 5) that result in k = 0.85. Both 
systems disqualify the HTA/ZE for its violating positivity. 
Unit of fcsT = 1 is used. 



To further demonstrate this property of the HQME, 
we report in Figs.|3]-[5] the results of transfer dynamics 
for three pairs of systems that are chosen rather arbitrar- 
ily. Two systems in a pair share a common (A, 7), and 
thus, are associated with one value of indicated in 
each figure. These results all indicate that the proposed 
HQME supports a broad range of applicability, which 
can be fairly well described by the quality parameter k 
of Eq. ([55|. 



All six cases presented in Figs. HHH are chosen to have 
the HTA/ZE unphysical, falling outside of its P-region, 
cf. Fig.[T] The corresponding HTA/ZE dynamics are 
therefore all wrong qualitatively and not shown. Rather 
we compare the results with an alternative HQME con- 
struction, labeled as M-HQME. This alternative is based 
on a similar approximation of the bath correlation func- 
tion, but treated on the Matsubara series expansion 
of Eq. (|XT|) . The resulting M-HQME differs from the 
present HQME only by the values of parameters Cr and 
A. Instead of Eq. dHj), the M-HQME assumes Cr = 
Re Co = A7C0t(^) and A = Aq = \[^ - cot(^)], 
from Eqs. (jA.2[) and (IA.5p . respectively. Apparently, the 
present HQME [Eq. (HUl) with Eq. (O], not just removes 
the singularity of cotangent function, but also performs 
better in dynamics, as compared to the exact results. 
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1.0F 




FIG. 5: Same as Fig.H but for {E° ,V) = (-2.5,3) (lower) 
and (E° , V) — (—8.5, 5) (upper), respectively, with same bath 
parameters of (A, 7) = (15,2) that result in k = 1.41 (same 
as the inset of Fig.[2ll. Both systems disqualify the HTA/ZE 
for its violating positivity. Unit of fcsT = 1 is used. 
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FIG. 6: Same as Fig.H but for {E°,V) = (-6,2.5) (lower) 
and {E°,V) = ( — 2,4) (upper), respectively, with same bath 
parameters of (A, 7) = (3, 4) that result in k = 2.38 (same 
as the main panel of Fig.[2|. Both systems disqualify the 
HTA/ZE for its violating positivity. Unit of ksT — 1 is used. 



B. Justification of the quality parameter 

Let us start with some insights for the qualitative su- 
periority of HQME over HTA/ZE. Note that the only 
approximation involved in each of them individually is 
the semiclassical bath correlation function Csc (i) [Eq. ^] 
versus its high-temperature version, Chta {t) [Eq. ([2])] , re- 
spectively. However, the resulting HQME acquires qual- 
itatively distinct modifications, especially with iC = 
iC + 5TZ appearing in Eq. PU)) , that lead to its being 
at least qualitatively consistent with the exact theory 
[Eq. (|A.6P ] . Physically, 5TZ serves as a diffusion correc- 
tion to the effective system Liouvillian, due to the ap- 
proximated treatment of bath correlation function. The 
HTA/ZE scheme, where STZ. = 0, does not have this 
diffusion modification to the system Liouvillian. Con- 
sequently, the applicability range of HTA/ZE depends 
not just on bath correlation function CmAit), but also 
sensitively on system. 

The HQME dynamics that includes explicitly the dif- 
fusion modified £' may thus support a system-insensitive 
criterion for its range of applicability. This argument has 
in fact been verified extensively in our numerical study 
of transfer systems. The applicability of HQME may 
therefore be addressed by examining the approximation 
involved in Csc(i). We quantify the approximation with 



the discrepancy function in the frequency-domain, 

SCiu) = Ciu) - [Csc{io) - A] 

= ^(cotl4--f + ^). (36) 
o;^ + 7"^ V 2 puj buj / 

The Drude model of Eq. ^ and A = /3A7/6 are used 
explicitly in the last identity. The discrepancy function is 
positive, symmetric, and monotonic decreasing, with the 
limiting values of 6C{uj = 0) = A and 6C{lu — > cxd) — 0. 
It is demonstrated in the inset of Fig. [71 Now let r(7) be 
the half-width at half-maximum, determined via 

<5C(c.)|^=r(,) - Y = (37) 

It can be well approximated by Eq. (j35bp . The result- 
ing r(7)/7, as participated in Eq. ((35)) . is demonstrated 
in Fig. [71 with no visible difference from the numerically 
exact evaluation of Eq. (|57|l . 

The criterion on applicability of HQME can now be 
considered for the condition under which SC{uj) can be 
treated as Markovian white noise for its effect of STZ 
on dissipative systems. We adopt the Kubo 's modu la- 
tion parameter-iSiiS for the criterion: k = y/T{j)/A = 
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FIG. 7: The scaled half- width-half-maximum, r(7)/7, as 
function of Pj. The approximate curve via Eq. (|35bp is found 
indistinguishable from the accurate one via Eq. p7|l . Shown 
in the inset are SC{lij)/A [cf. Eq. (|36p ] versus f3ij, at some se- 
lected values of 7 (unit of ksT). The curve does not change 
for /37 < 1. 
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FIG. 8: Validity diagram of HQME by the criterion of the 
parameter k [Eq. (|35p ] . 



a = 33 (the lower panel of Fig. [3]). 



VI. CONCLUDING REMARKS 

We have proposed the HQME that may be used with 
good confidence to arbitrary systems. The theoretical de- 
velopment is rooted in an improved semiclassical treat- 
ment of Drude bath. This alone improves the conven- 
tional high-temperature or classical approximation by 
two or three orders in temperature parameter, respec- 
tively, as argued in SecUl The resulting HQME can 
be considered as a natural extension and modification of 
the conventional stochastic Liouville equation and Zus- 
man equation. While it retains their appealing physi- 
cal pictures and numerical efficiency, the HQME shows 
remarkable improvement over the conventional theories. 
Its broad range of validity and applicability, in terms 
of density matrix positivity and dynamics quality, are 
extensively demonstrated on two-level model systems. 
We have also proposed a criterion to estimate the per- 
formance of HQME. This criterion depends only on the 
system-bath coupling strength, characteristic bath mem- 
ory time and temperature. Our results all reveal that the 
HQME may serve as a versatile tool wherever the exact 
approach is numerically too expensive to afford. 
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•\/6r(7)/(/3A7) > 1. This is the quality parameter of 
Eq. dSSD. 

The range of applicability of HQME, as depicted by the 
region of k > 1 in Fig.[8l is quite impressive. Consider, 
for example, the fact that it covers the value of /3A < 57, 
when /37 < 1. For a rather fast solvation of = 100 fs, 
that /?7 < 1 would support temperature T > 75 K and re- 
organization energy at least the range of A < 3000 cm^^. 
This covers almost all systems of practical interest. Note 
that all dynamics results presented in this work are in 
the strong system-bath coupling regime. Let a be the 
dimensionless system-bath coupling strength parameter, 
defined via [J(w)/w](^=o = 7ra/2, as the Kondo param- 
eter in the Ohmic model. We have a = 4A/(7r7). The 
present HQME would support a < 72, with 7"^ = 100 
fs and T = 75 K, if K > 1 criterion is used. The cal- 
culations demonstrated in this paper have the coupling 
strength ranging from a = 0.95 (Fig. [2] and Fig. (5]) to 



APPENDIX: EXACT HEOM FORMALISM 

The HEOM formalism for the reduced system dy- 
namics can be constructed using the calculus-on-path- 
integral technique i^i^ili^i^ i ^ 1,12^13 ^pj^jg formalism is in 
principle exact and nonperturbative. However, its spec- 
ified form depends on the way of treating the bath 
correlation function, under the constraint of the exact 
fluctuation-dissipation theorem [Eq. (|3])]. Consider the 
Matsubara expansion method, in which the bath corre- 
lation function with the Drude model of Eq. ([3]) reads 

00 

C(t>0)=coe-^* + ^Cfce-'^'=*. (A.l) 

The flrst term, with 

co = A7[cot(/37/2)-i] , (A.2) 
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arises from the Drude pole. The sum term arises from the 
Matsubara poles or frequencies of 7fe>i = 2fc7r//3, with 



8fc7rA7 



(2/c7r)2 - (/37)2 



fc > 1. 



(A.3) 



To construct the HEOM, the infinite sum over the k 
index in Eq. (jA.ip need to be truncated. To that end, let 
7o = 7 and recast Eq. (jA.ip by 



(A.4) 



k=0 



where 



Cfc 



k=K+l "^'^ 



/37 ^ 2 \ 



K 



This treatment is in principle exact if the K is chosen 
large enough and the resulting reduced system density 
matrix dynamics of primary interest is converged. 
The resulting HEOM formalism reads^iiiiia^ 



k=Q 
K 

i E v/("fe + l)|cfe| [Q, 

fe=0 
K 

i E V"fe/|cfe|(cfeQp, 



(A.6) 



fc=0 



with 



STZKO^AKiQdQ.O]]. 



(A.7) 



The reduced density operator of primary interest is 
p = po- The subscript n = {rik > 0; k — 0,---,K}, 
which consists of a set of nonnegative indices, specifics 
in general a given auxiliary density operator (ADO) of 
Pn = pno---nK- The subscript differs from n only by 
changing the specified to rife ± 1. An iV"^-tier ADO is 
of no + • • • + riK = N, and the total number of ADOs at 
this tier is ^^t^?' ■ Each p„ in Eq. (|A.6p has been scaled 
individually to have a uniform dimension and error toler- 
ance. This validates a simple, on-the-fly, filtering algo- 
rithm for efhcient numerical propagation, including au- 
tomatic truncation of hierarchy.^^'^^ However, the exact 
HEOM calculations is still numerically expensive, espe- 
cially for complex systems, under strong coupling with 
bath at low temperature. 
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